Exact solution of stochastic directed sandpile model 
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We introduce and analytically solve a directed sandpile model with stochastic toppling rules. The 
model clearly belongs to a different universality class from its counterpart with deterministic toppling 
rules, previously solved by Dhar and Ramaswamy. The critical exponents are Du = 7/4, r = 10/7 
in two dimensions and Dii = 3/2, r = 4/3 in one dimension. The upper critical dimension of the 
model is three, at which the exponents apart from logarithmic corrections reach their mean-field 
values Dm — 2, t — 3/2. 
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Numerical and analytical studies of sandpile models 
of Self-Organized Criticality [Q continue to be a subject 
of considerable research activity. In particular, a lot of 
effort has been recently invested in establishing the set of 
universality classes in these systems 0. The consensus 
seems to be that the universality class of a d-dimensional 
sandpile model is determined by answers to the following 
list of questions: 

(i) Is it a critical slope or a critical height model? 
In other words, does a site topple when its local slope 
or its height exceeds a certain threshold value. Critical 
height models (e.g. Bak-Tang-Wiesenfeld (BTW) model 
[jlj ) were studied more extensively in the past and are in 
general better understood. 

(ii) Is sand redistributed isotropically in a toppling 
event? According to this property sandpiles could be 
classified as isotropic or directed (anisotropic) . The com- 
mon knowledge is that it is a relevant parameter, i.e. an 
arbitrary small anisotropy of toppling rules usually drives 
the model to the directed universality class. 

(iii) Finally, it is important if the sand is distributed 
deterministically or randomly in each individual top- 
pling. In a model with deterministic toppling rules the 
configuration of the bulk of the sandpile remains un- 
changed if every single site on the lattice topples exactly 
once. This additional symmetry is usually important for 
the universality class of the model: e.g. the determinis- 
tic one-dimensional isotropic critical height sandpile (ID 
BTW) has only trivially distributed avalanches of frac- 
tal dimension 2, while its cousins with additional ran- 
domness in toppling rules, such as the Zaitsev model ||, 
Oslo model Linear Interface Model Q, etc., seem to 
belong to the same universality class where avalanches 
have a non-integer fractal dimension D ~ 2.23 and a 
probability distribution with a clean power law exponent 
t ~ 1.27. 

Despite much numerical and analytical effort on the 
original BTW sandpile model (which 
is a deterministic/isotropic/critical-height model in the 
above classification), its critical exponents in two dimen- 



sions still remain controversial 0] . The situation is some- 
what better for directed models. Soon after the original 
BTW sandpile model Dhar and Ramaswamy intro- 
duced and exactly solved in all dimensions its directed 
cousin - the Dhar-Ramaswamy (DR) model [|| . 

Both BTW and DR models have deterministic toppling 
rules. As far as stochastic models are concerned there is 
preciously little analytical results. Apart from an ex- 
act solution of a model, equivalent to the ID stochastic 
directed sandpile |t]], stochastic sandpiles were studied 
only numerically. In this paper we present an analyti- 
cal study of the stochastic directed sandpile model in all 
dimensions. Stochastic directed sandpiles were brought 
to the attention of the community in two recent papers 
, reporting numerical studies of several variants of such 
models in two dimensions. During the preparation of this 
manuscript, there appeared a closely related preprint by 
Paczuski and Bassler in which an analytical study of 
the directed stochastic sandpile model was presented and 
similar results were obtained. In particular, using differ- 
ent analytical arguments they have arrived at the same 
set of exponents. 

The microscopic rules of the stochastic directed sand- 
pile model that we selected to study are closely related 
to those of the DR model || . These rules are modified in 
the spirit of a stochastic isotropic sandpile model known 
as the Manna model jl(j. It is easier to define our rules 
in two dimensions, while generalization to higher dimen- 
sions is straightforward. A stable configuration of our 
model is specified by the integer height of the sandpile 
z{xx,x<z) < 1 at each point of a 2D square lattice. The 
lattice has open boundary conditions along the diagonal 
coordinate, x\\ = x\ + x%, and periodic boundary con- 
ditions in the transversal direction x± = x± — x%. The 
sand is added randomly at the line with x\\ = and falls 
off the edge at zii = Ln. The difference between our 
model and the DR model lies in toppling rules. In both 
cases, once the height at any given site exceeds one, this 
site becomes unstable and loses two grains of sand to its 
nearest neighbors in the direction of increasing xm. How- 
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ever, while in the DR model each of these two neighbors 
gets exactly one grain of sand, in our stochastic variant 
the decision where to move any particular grain is done 
independently for each grain. In other words, with prob- 
ability 1/4 both grains end up on the left neighbor, with 
probability 1/4 they go to the right neighbor, and only 
with probability 1/2 will each neighbor get one grain as 
in the DR model. Obviously, on average each neighbor 
gets one grain, yet the additional stochastic element in 
these rules drives the model away from the universality 
class of the DR model . It is easy to see that unlike the 
deterministic rules of the DR model, the new stochastic 
rules allow for multiple topplings of some sites within one 
avalanche. Indeed, let us consider an example where the 
first active site of the avalanche toppled one grain of sand 
to each of its two neighbors. Let us further assume that 
these neighbors both had z = 1 so that they both top- 
pled and by chance distributed all four of the resulting 
grains of sand to the same site in the next layer. This 
site has received four grains of sand and is guaranteed 
to topple twice. The numerical simulations (|] confirm 
the existence of multiple topplings in other variants of a 
stochastic directed sandpile model. 

In order to get an analytical handle on the properties 
of our model we employ the same trick that was used 
by one of us to solve the ID directed stochastic sandpile 
model j7|. Due to the Abelian nature of the model, we 
can change the order in which topplings are performed 
without changing the outcome. It is convenient to do 
topplings layer by layer. This means that we topple any 
given site as many times as necessary to make it stable 
before moving on to the next unstable site, and we topple 
all unstable sites in one layer (a given xn) before toppling 
any site in the next layer. Let us concentrate on a site 
with coordinates xn and x± immediately after we have 
finished with topplings in the (xii — l)-th layer. Assume 
that two of its neighbors with coordinates x\\ — 1 and 
x±_ ± 1 have toppled, respectively, n\ = n(xu — 1, x± — 1) 
and ni = n(x\\ — 1, x± + 1) times. The average number of 
grains of sand that our selected site would receive from 
the previous layer is (n\ + 712). In the DR model there 
are no fluctuations around this average. Also, due to 
the absence of multiple topplings in this deterministic 
directed model, ri\ and ni can be only or 1. Therefore, 
in the DR model a site can receive either n\ + ni = 2 
grains, in which case it is guaranteed to topple exactly 
once, or n\ + ni = 1, in which case it can topple with 
probability 1/2 (i.e. it topples if it had z — 1 before the 
transfer and remains stable if it had z = 0). From this 
one can show Q that in the DR model the set of sites 
which topple at each layer form an interval with no holes 
inside. The size of this interval as a function of the layer 
number x\\ performs an ordinary random walk. 

In the stochastic model the relation between the num- 
ber of topplings in two subsequent layers is more com- 
plicated. Let us concentrate on the behavior of the total 



number of topplings N(x\\) = n { x \\t x ±) m a given 
layer xii. This number fixes the number of grains of sand 
transferred from the layer xii to the next layer as 2iV(x||). 
It is easy to see that a site which has received an even 
number 2k of grains of sand from a previous layer will al- 
ways topple exactly k times and, therefore, will transfer 
the same 2k grains of sand to the layer directly below it. 
That means that as far as N(x\\) is concerned, such sites 
behave in a completely passive manner, i.e. they do not 
lead to a decrease or an increase of the total number of 
topplings N(x\\) from layer to layer. On the other hand, 
any site which received an odd number 2k + 1 of grains of 
sand from a previous layer has equal chance to topple k 
times (if it had z — before the transfer) or k+1 times (if 
it had z = 1). In the former case this site would decrease 
the grain flow 2 A^(a; 1 1 ) by one, while in the latter - in- 
crease by 1. Let us call any site which has received an odd 
number of grains of sand from the previous layer an ac- 
tive site. The equation, which is the central result of this 
work, relating the change in the total number of topplings 
from layer to layer to the number of active sites AT a (x||) in 

a given layer is N(xu) — N(xu — l) + i So^i* Co; where 
all £ a are —1 or +1 with equal probability and indepen- 
dent of each other. These random numbers correspond 
to whether each of the N a (x\\) active sites had the height 
z = or z = 1 before the avalanche started. It is straight- 
forward to demonstrate that, as in the DR model, in the 
steady state of the directed stochastic model all possible 
stable configurations of z are equally represented, and, 
therefore, there are no correlations between the heights 
at different sites, and each height is equally likely to have 
z = or z = 1. It is more convenient to rewrite the above 
equation in a continuous notation, which works as long 
as iV Q (x||) >> 1 : 



dN(xu) 1 / : — - , , 



(1) 



Here n(t) is a standard Gaussian variable with zero mean 
and the standard deviation equal to unity. This equation 
describes an unbiased random walk iV(x||) vs xi with a 
variable step, given by ^yjN a (x\\). A random walk (an 
avalanche) starts with A^(0) = 1 and ends at x\\ when 
N(x\\) < for the first time. Let us assume that both 
N(x\\) and N a (x\\) in a surviving avalanche scale with X|| 
with the exponents a and a a , respectively. Alternatively 
one can say that they are related to each other by N a ~ 
N aa / a . Plugging this relation into Eq. ([!]), after easy 
algebra we get the exponent relation 



a 



(2) 



In addition to this exponent relation, the mapping of 
the avalanche process onto a random walk immediately 
gives us the power law T|| in the probability distribu- 
tion of avalanche sizes. Indeed, an avalanche ends when 
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the random walk described by Eq. (Q) first enters the 
N(x\\) < semi-axis, and it is a standard result of the 
theory of stochastic processes that the distribution of 
first returns of a generalized random walk has an ex- 
ponent T|| = 1 + a. This agrees with the well-known 
exponent relation th = Dm valid for a general directed 
sandpile model. Indeed, the fractal dimension Dm, de- 
fined by Yfi=i N(i) = S ~ x, " , is obviously related to a 
through Dm = 1 + a. 

The Eq. ([!]) applies equally well to the DR and the 
stochastic directed sandpile models. The difference be- 
tween these two models lies only in the scaling of the 
number of active sites with xii. As was explained above, 
in the 2D DR model the only two active sites lie at the 
edge of an interval of toppled sites. Indeed, only these 
sites get 1 grain of sand, while the rest get either or 
2. Therefore, in the 2D DR model iV a (x||) = 2 is just 
a constant, a a = 0, and the Eq. (fy) describes an ordi- 
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nary random walk, in which N(x\\) ~ xjj = in . The 
introduction of a stochastic element in particle redistri- 
bution dramatically changes the number of active sites 
at any given layer of the avalanche. Indeed, when grains 
arc distributed independently, any site which has at least 
one toppled neighbor in the previous layer is equally likely 
to receive an even or odd number of grains of sand, 
and therefore, it has a probability 1/2 of becoming ac- 
tive. Thus, in the stochastic model the exponent a a 
defines how the number of distinct sites that topple at 
least once, scales with the layer number ajit. The dif- 
ference between exponents a and a a comes solely from 
the existence of multiple topplings. These two exponents 
have to obey the inequality a > a a , and their difference 
a — a a determines how the average number of topplings 
ntop{x\\) at a given site in Xii-th layer scales with X||: 
n top ~N/(2N a ) ~a;fp a «. 

We proceed with an argument that in the 2D directed 
stochastic model a a = 1/2, and, therefore, by the virtue 
of Eq. (Q) a = 3/4. It is a straightforward task to 
determine the average number of topplings (n(x\\, x±)) 
at a given site xm,x_l, where the average is performed 
over the whole ensemble of avalanches, so that avalanches 
that die out before reaching this site, contribute to the 
average. As was noted in due to the conservation of 
sand and the stationarity of the sandpile (n(xn, x±)) has 
to satisfy the diffusion equation with a source: 

d(n(x\\,x ± )) 1 d 2 {n{x ]bXl _)) 

^ =2 dxl W 

This average balance equation is also exact for our 
stochastic model, where it proves that, like in the DR 
model, sites that topple at least once are spread over 
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the interval of length Ax± ~ xJ in the X|| layer. In 
the DR model these sites form a dense interval with no 
holes, and, therefore, their number is known to scale 



exactly as xJ . The situation is somewhat less obvi- 
ous in the stochastic directed model, where the set of 
toppled sites can have holes. However, one can argue 
that these holes would mostly be concentrated near the 
boundaries of the avalanche in any given layer, while the 
majority of toppled sites in the center would form a com- 
pact interval. Indeed, as will be confirmed later, the 
2D stochastic directed model is characterized by multiple 
topplings, where a site at a layer x\\ would typically top- 
ple ntop(x||) ~ x^ 4 times within one avalanche. Since 
any of the 2ntop grains is equally likely to go to each of 
the two nearest neighbors in the next layer, the situa- 
tion where one of these neighbors would receive less than 
two grains is exponentially unlikely for large n top . But 
once a site has received two or more grains of sand it is 
guaranteed to topple at least once. Therefore, both near- 
est neighbors down the slope from the site which toppled 
many times would most likely topple at least once. In 
other words, the creation of a new hole (a region free of 
topplings) is exponentially suppressed near sites which 
themselves toppled many times. This means that for a 
sufficiently large xi the majority of sites (especially those 
close to the center of the avalanche region) will belong to 
a hole-free region. Since the size of the interval covered 
by an avalanche scales as x^ 2 , the number of active sites 
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should also scale as N a ~ x,, . From a a = 1/2 with the 
help of Eq. (Q) one gets a = 3/4, t\\ = D\\ = 7/4, and 
t = 1 + (t|| — 1)/Dii = 10/7. These results are in a nice 
agreement both with previous numerical simulations of 
various versions of stochastic directed sandpile model in 
two dimensions B and with our own simulations of the 
model. In Fig. 1 we present the results of our simula- 
tions for the effective exponents a = dlogN/dlogx^ and 
a a = d\og N a / d\ogx\\ as a function of xii. The numeri- 
cal exponent a nicely agrees with the analytical results. 
The exponent a a is less clean due to the presence of holes 
near the boundary of the avalanche region. The exponent 
seems first to overshoot to a value of almost 0.6 but then 
goes down so that in the end of the range of our simu- 
lations, xi ~ 30000, it is consistent with our theoretical 
prediction a a = 1/2. 

Unlike its deterministic cousin, the stochastic directed 
sandpile model exhibits a non-trivial scaling even in one 
dimension. The ID deterministic directed sandpile model 
is trivial in the sense that there is just one SOC configura- 
tion and the addition of a grain of sand always results in 
an avalanche of L\\ topplings in which this grain is trans- 
ported and discarded at the open boundary of the system. 
A ID stochastic directed sandpile model can have several 
variants of simple microscopic toppling rules. In one vari- 
ant, which is essentially identical to the model studied by 
one of us in |7j] , once a height at a given site exceeds one, 
either one or two grains are transferred to the nearest 
neighbor down the slope. It is easy to see (for details see 
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) that this model is equivalent to a ID random walk so 
that N a = const, while the typical number of topplings 
N scales as a function of Xii — x as N(x) ~ x 1 / 2 . The 
distribution of avalanche spatial length in this model has 
an exponent t\\ — 3/2, while that of avalanche volume - 
t = 4/3. In another variant of the ID stochastic toppling 
rules an unstable site always loses two grains of sand and 
each of these grains with equal probability goes to the 
nearest or next-nearest neighbors down the slope. 

As in the DR model the upper critical dimension for 
the stochastic directed sandpile model is d u = 3. In 
this dimension the expected number of topplings at each 
site in a layer x\\ grows only logarithmically with x\\. 
Therefore, a — a a apart from the logarithmic correc- 
tions. From Eq. (0) in this case we get a = (1 + a)/2, 
which has the solution a = 1, rn = 2, r = 3/2. This is 
a standard set of mean-field exponents for any branch- 
ing (avalanche) process in high enough dimension. In 
Fig. 2 we plot the numerical effective exponents in the 
3D stochastic directed model. They agree well with the 
mean field values. 

In conclusion, we have found an analytic solution of 
the stochastic directed sandpile model in any dimension. 
The main difference of this model from its deterministic 
counterpart - the Dhar-Ramaswamy model - lies in the 
fractal dimension of the set of active sites, i.e sites that 
can add or remove one grain from the overall flow of sand 
between two subsequent layers. Whereas in the 2D DR 
model in any layer there are only two sites at the edges 
of the interval of toppled sites which are active, in the 2D 
stochastic directed sandpile model each of approximately 
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af|l toppled sites in this interval has a 1/2 chance of be- 
ing active in the above sense. This leads to an increase in 
the fractal dimension of an avalanche from D\\ = 3/2 to 
Du = 7/4 due to the appearance of multiple topplings. 
The difference between critical properties of stochastic 
and deterministic directed models disappears in high di- 
mensions d > 3, where multiple topplings in a stochastic 
directed sandpile become prohibitively unlikely and all 
exponents acquire their mean-field values. 
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FIG. 1. The scale dependent effective exponents 

a = dlog N / d\ogx\\ and a a — dlog iV a /dlog X| as a func- 
tion of xu for the 2D model, and for two longitudinal system 
sizes L|| = 5000 and L\\ = 40000. 
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FIG. 2. (a) The effective exponents D,a,a a , and r as a 
function of x\\ for the 3D model, (b) The expected number 
of topplings at each site n t0 p as a function of x\\. Note the 
logarithmical dependence on xu. 
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